%% Step 3

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
% Author:Produced for JdG,AL,IM by Christopher Evans at UPF 
%
% Program: Market clearing condition to calculate the wedges in the model. 
%          Once wedges have been constructed then we will loop over the
%          whole model to calculate X_hat (expenditure hat) (step 4),
%          the change needed for the wedges to be zero. Once we have the 
%          value of X_hat with no shocks in the model we can use it for 
%          future counterfactual exercises.
%
%
% Differences with expenditure function: Step 3 calculates the wedges. The
% expenditure function loops over a market clearing condition with the
% wedge set to 0 until the model converges at a new equilibrium. 
%
% NEW: For all code, we will change the subscripts to match those in
% master_notes3.pdf and the draft. I.e., we will change to {mn,ij} pairs,
% which represent country pair {mn} and sector pair {ij}. Global variables 
% changed and slight modification to file structure for calling up data.
% We also get rid of gamma_type variable. What previously _2.csv. files
% slightly as well as adjusting in our new directory.
%
% We also get rid of using the .(ISO{}) cell call up for the gammas and
% instead use a large matrix saved in Step2. This will help save on loops
% and speed up the procedure.
%
% Last Updated: 5/02/2019 by IM & JDG
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clearvars

global alphaT rhoT lambdaT etaT rho sigma lambda eta varphi ...
    tol maxit Pi_hat_n D_hat_n Xi_hat_mnjf Xi_hat_mnj a_hat_f a_hat ...
    N J FI_J france countrysecD firmsecD_sorted start start_sorted ...
    sum_by_country_dummy ISO ... 
    sL_n sPi_n sD_n alpha_WIOD_j alpha_WIOD_francej labour_share_PWT
 

% Loading the MAT files of step 1 and step 2:

load MAT/WIOD.mat
eval([strcat('load MAT/Step_2_firm_data_alpha',num2str(alphaT),'.mat')]);


%% All countries except france, sector level

% X_pi: Need to multiply by  X pi here and sum over the countries n: 

X_pi= sum(pi_mnj.*X_mnj,2);

% Making the previous result diagonal so that we can multiply it by the
% gamma matrix. Diagonal so elements of value on the diagonal and the rest
% are zero

diag_X_pi=diag(X_pi);

for m=1:N    
   for n=1:N
      %for a given i and n gamma is a [JxJ] matrix, the X_pi is [J x1]
      %but if we diagonalize it will be a [JXJ] wilth the value along
      %the diagonal and zeros elsewhere. 
        
      gamma_X_pi(start(m):start(m+1)-1,start(n):start(n+1)-1)= ...
         gamma(start(m):start(m+1)-1,start(n):start(n+1)-1)*...
         diag_X_pi(start(n):start(n+1)-1,start(n):start(n+1)-1);   
    end
end


if alphaT == 1 | alphaT == 2
   % oneminus_alpha: 1-labor share, which needs to be scaled  up to JN*JN
   % levels. This is done by using the countrysec dummy and then repmat to
   % extend the number of columns.

   oneminus_alpha_vector=countrysecD*(1-labour_share_PWT);

   oneminus_alpha_matrix=repmat(oneminus_alpha_vector,[1 N*J]);
else
   % oneminus_alpha: This one takes the alpha from the WIOD; by sector and
   %                scales up to the JN*JN level 

   oneminus_alpha_vector=(1-alpha_WIOD_nj)';
   oneminus_alpha_matrix=repmat(oneminus_alpha_vector,[1 N*J]);
end

% one_minus_alpha_gamma_X_pi: Multiply gamma by 1-alpha

one_minus_alpha_gamma_X_pi= gamma_X_pi.*oneminus_alpha_matrix';

% rho_nj_repmat: Expanding rho by importing countries n 

rho_nj_repmat=repmat(rho,[N 1]);

% rho_mnij_repmat: Expanding rho to fit the WIOD matrix [N*J x N*J]

rho_mnjj_repmat=repmat(rho_nj_repmat',[N*J 1]);

rho_one_minus_alpha_gamma_X_pi=((rho_mnjj_repmat-1)./rho_mnjj_repmat).*...
    one_minus_alpha_gamma_X_pi;

  
                        
% Summing over i, meaning we should have a [NJ N] afterwards 

for i=1:N
   %loops through each country m summing the i sectors 
   
   pseudo_intermediates(:,i)= sum(rho_one_minus_alpha_gamma_X_pi(:,start(i):start(i+1)-1), 2);
end

% Wedge: Difference in actual exports from WIOD and pseudo exports

wedge =X_mnj-(pseudo_intermediates+PC_mnj);

X_mnj=X_mnj+0.0000000000000001.*(X_mnj==0);

% Wedge over X 

wedge_over_X_mnj = wedge./X_mnj;

%% Wedge using French firms 

% X_mnjf: Multiplies the expenditure matrix for france with the firm dummy
%         to create a [F x N matrix]

X_mnjf = firmsecD_sorted*X_mnj(start(france):start(france+1)-1,:);

% pi_X_mnjf: Multiplying firm shares of expenditures (pi_mnj) by the
%            expenditure matrix that has been scaled up by the firm sector
%            dummy

pi_X_mnjf=pi_mnj_sorted.*X_mnjf;

% mu_pi_X_mnjf: Multiplying by inverse of the markup to recover the share
%               that pays inputs
repmat_sigma = repmat(firmsecD_sorted * sigma,[1,N]);
repmat_rho = repmat(firmsecD_sorted * rho,[1,N]);

mu_mnjf = repmat_sigma.*repmat_rho./(repmat_sigma.*(repmat_rho-1)-(repmat_rho-repmat_sigma).*pi_mnj_sorted);
mu_constant = repmat_rho./(repmat_rho-1);
plot(mu_constant,mu_mnjf);
mu_pi_X_mnjf = pi_X_mnjf ./ mu_mnjf;

%pi_X_njf: This is the sum over country k, should result in a [Fx1] vector

mu_pi_X_njf=sum(mu_pi_X_mnjf,2);

%one_minus_alpha_f: Labor shares for each firm = 1-alpha

one_minus_alpha_f=1-alpha_f_sorted;

% Multiplying the previous summed equation by this labor share. This is
% done before multiplying by gamma, which will come later. 

mu_pi_X_laborshare_njf=mu_pi_X_njf.*one_minus_alpha_f;


% mu_pi_X_laborshare_njf: Repurposing the matrix such that it spans across the
%                      j sectors of country m. It will have the same value
%                      j and make the following computation easier

mu_pi_X_laborshare_njf_repmat= repmat(mu_pi_X_laborshare_njf, [1 J*N]);     

mu_gamma_pi_X_laborshare_njf = mu_pi_X_laborshare_njf_repmat.*GAMMA_sorted_matrix;
french_Z_mnj = sum(mu_gamma_pi_X_laborshare_njf,1);

% stacked_french_Z_mnj: Stacking the previous result puts it in the same
%                       format as the rest of the Matlab file, such that we 
%                       will have a [J*N 1] vector, where the 1 is due to
%                       looking at only one n, france.

stacked_french_Z_mnj=reshape(french_Z_mnj, [N*J 1]);


% X_mnj: This is previously computed by is done again to avoid any issues
%        with previously adding epsilon (0.0000...01), which was done to 
%        avoid nan problems

X_mnj=Z_mnj+PC_mnj;

% X_in: Aggregating the expenditure matrix to the country level

for m=1:N
    X_mn(m,:)=sum(X_mnj(start(m):start(m+1)-1,:),1);
end

% french_PC_mnj: Final consumption when n=france

french_PC_mnj=PC_mnj(:,france);

% french_X_mnj: Exports from country m sector j to france

french_X_mnj=X_mnj(:,france);

% french_wedge_mnj: Calculating the wedge for france only

french_wedge_mnj=french_X_mnj-(stacked_french_Z_mnj+french_PC_mnj);

% Adding epsilon to reduce NAN issues 

french_X_mnj=french_X_mnj+0.0000000000000001.*(french_X_mnj==0);

% Calculating the wedge over X for france

french_wedge_over_X_mnj=french_wedge_mnj./french_X_mnj;

% Combining all countries wedge and France 

% zeta_mnj: Is the definition for the wedge

zeta_mnj=wedge;

% Replacing the wedge for n=france

zeta_mnj(:,france)=french_wedge_mnj;